Tecniche di data mining per la classificazione di soggetti fumatori

La presente analisi si propone di affrontare un problema di classificazione che ha come obiettivo discriminare tra individui fumatori e non sulla base di una serie di misurazioni mediche. Le tecniche di data mining applicate sono: alberi decisionali, KNN, regressione logistica e random forest.

I dati utilizzati possono essere scaricati al seguente link: https://www.kaggle.com/datasets/kukuroo3/body-signal-of-smoking

1 Introduzione

Il fumo è una delle più grandi minacce per la salute in quanto rappresenta uno dei maggiori fattori di rischio nello sviluppo di patologie neoplastiche, cardiovascolari e respiratorie. Secondo i dati diffusi dall’OMS è responsabile del decesso di almeno 6 milioni di persone ogni anno. Questo si traduce in ingenti spese sostenute dai servizi sanitari nazionali per le cure di pazienti affetti da patologie causate dal fumo.

L’esempio qui considerato è quello del National Health Insurance Service (NHIS), sistema di assicurazione sanitaria nazionale del Sud Corea. Nel 2014 il NHIS ha indetto una causa contro KT&G, il principale produttore di tabacco coreano, richiedendo un risarcimento di 53,3 miliardi di won, pari alla spesa totale sostenuta dal 2003 al 2013 per il trattamento di pazienti affetti da cancro ai polmoni e alla laringe, fortemente legati al fumo.

L’NHIS ha raccolto varie informazioni sui check-up svolti da un campione di cittadini assicurati dal servizio stesso. Queste sono costituite da dati identificativi dei pazienti, dai risultati della visita e dallo stato di fumatore.

La seguente analisi applicherà quindi diversi metodi di classificazione al fine di identificare lo stato di fumatore per pazienti di cui si conoscono solo valori biologici. Questo permetterebbe all’NHIS di avere una previsione su quanti tra i loro assicurati siano fumatori o meno.

2 Dataset

Il dataset contiene 55692 osservazioni. La variabile target Smoking assume valore 0 nel caso di individuo non fumatore e 1 in caso di individuo fumatore. Il 63% degli individui è non fumatore, il 37% lo è. Sono inoltre presenti 26 variabili così definite:

1. ID: codice identificativo.
2. Gender: genere.
3. Age: età (anni); rilevata arrotondando a valori multipli di 5.
4. Heigh: altezza (cm); misurata arrotondando a valori multipli di 5.
5. Weight: peso (kg); misurata arrotondando a valori multipli di 5.
6. Waist: circonferenza della vita (cm).
7. Eyesight (left): vista occhio sinistro.
8. Eyesight (right): vista occhio destro.
9. Hearing (left): udito normale (0) o presenza di anomalie (1) nell’orrecchio sinistro.
10. Hearing (right): udito normale (0) o presenza di anomalie (1) nell’orrecchio destro.
11. Systolic: pressione sistolica (massima), misurata durante la contrazione ventricolare.
12. Relaxation: pressione diastolica (minima), misurata durante il rilassamento ventricolare.
13. Fasting blood sugar: glicemia, ossia il valore di glucosio nel sangue.
14. Cholesterol: colesterolo totale, molecola lipidica in parte sintetizzato nel fegato e in parte introdotto con l’alimentazione.
15. Triglyceride: trigliceridi, lipidi introdotti principalmente con la dieta marcatori della salute metabolica.
16. HDL: lipoproteine che ripuliscono le arterie dal colesterolo in eccesso trasportandolo al fegato, dove viene smaltito (“colesterolo buono”).
17. LDL: lipoproteine che trasportano il colesterolo dal fegato alle cellule del corpo e che, in quantità eccessiva, si depositano sulle pareti arteriose provocandone un progressivo inspessimento (“colesterolo cattivo”).
18. Hemoglobin: emoglobina, proteina contenuta nei globuli rossi deputata al trasporto di ossigeno nel sangue.
19. Urine protein: proteine nelle urine, valori elevati possono indicare una malattia renale; assume valori tra 1 e 6.
20. Serum creatinine: creatinina, sostanza di scarto prodotta nei muscoli indice della salute dei reni.
21. AST: aspartato transaminasi, enzima il cui aumento del livello può indicare una patologia o lesione al fegato.
22. ALT: alanina amino transferasi, enzima epatico che aiuta l’organismo a metabolizzare le proteine, anch’esso indicatore della salute del fegato.
23. γ-GTP: gamma glutamil transferasi, enzima epatico le cui concentrazioni elevate possono indicare una lesione al fegato o ai dotti biliari.
24. Oral: indica se il paziente ha svolto una visita dentistica (1) o meno (0).
25. Dental caries: presenza (1) o assenza (0) di carie.
26. Tartar: presenza (yes) o assenza (no) di tartaro.

smoke <- read.csv("smoking.csv")
smoke$gender<-as.factor(smoke$gender)
smoke$hearing.left.<-as.factor(smoke$hearing.left.)
smoke$hearing.right.<-as.factor(smoke$hearing.right.)
smoke$dental.caries<-as.factor(smoke$dental.caries)
smoke$tartar<-as.factor(smoke$tartar)
smoke$smoking<-as.factor(smoke$smoking)

Partendo dal dataset originale appare opportuno procedere subito all’eliminazione di alcune variabili, quali

• ID: inutile ai fini dell’analisi.
• Oral: assume valore 1 per tutte le osservazioni e quindi non aggiunge informazione.
• Urine protein: la fonte non fornisce informazioni sui 6 livelli assumibili dalla variabile ed essi non sono quindi riconducibili al valore numerico della concentrazione di proteine nelle urine. Per non inficiare le analisi è stato scelto di eliminare la variabile, avendo precedentemente considerato che il set di variabili rimanenti è opportunamente in grado di spiegare la variabile target.
• Eyesight left e eyesight right: la fonte non fornisce informazioni precise sul metodo di rilevazione della vista e sul significato dei valori rilevati. Conducendo un’analisi più specifica, in particolare concentrandosi sui boxplot, si evidenzia una distribuzione anomala: le osservazioni sono concentrate tra 0.1 e 2 e appaiono poi molti valori pari a 9.9, lasciando un gap illogico. Il fatto che ci sia una polarizzazione su un numero così fuori dal range non è in alcun modo spiegabile dalle informazioni in nostro possesso, quindi verranno rimosse le variabili per non inficiare le analisi successive. Inoltre, i boxplot condizionati alla classe sono identici, quindi presumibilmente non vanno a discriminare sullo stato di fumatore.

ok <- ggplotly(ggplot(smoke, aes_string(x = "smoking", y = smoke$eyesight.right., col = "smoking", fill = "smoking"))+geom_boxplot(alpha=0.2)+ylab("eyesight.right"))

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

ok2 <- ggplotly(ggplot(smoke, aes_string(x = "smoking", y = smoke$eyesight.left., col = "smoking", fill = "smoking"))+geom_boxplot(alpha=0.2)+ylab("eyesight.left")+labs(title="boxplot di eyesight.right sopra e di eyesight.left sotto"))
subplot(ok,ok2,nrows=2)

library(dplyr)
smoke<-smoke%>%dplyr::select(-c(ID,eyesight.left., eyesight.right., Urine.protein, oral)) #rimozioni variabili 

str(smoke)

## 'data.frame':    55692 obs. of  22 variables:
##  $ gender             : Factor w/ 2 levels "F","M": 1 1 2 2 1 2 2 2 1 2 ...
##  $ age                : int  40 40 55 40 40 30 40 45 50 45 ...
##  $ height.cm.         : int  155 160 170 165 155 180 160 165 150 175 ...
##  $ weight.kg.         : int  60 60 60 70 60 75 60 90 60 75 ...
##  $ waist.cm.          : num  81.3 81 80 88 86 85 85.5 96 85 89 ...
##  $ hearing.left.      : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ hearing.right.     : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ systolic           : num  114 119 138 100 120 128 116 153 115 113 ...
##  $ relaxation         : num  73 70 86 60 74 76 82 96 74 64 ...
##  $ fasting.blood.sugar: num  94 130 89 96 80 95 94 158 86 94 ...
##  $ Cholesterol        : num  215 192 242 322 184 217 226 222 210 198 ...
##  $ triglyceride       : num  82 115 182 254 74 199 68 269 66 147 ...
##  $ HDL                : num  73 42 55 45 62 48 55 34 48 43 ...
##  $ LDL                : num  126 127 151 226 107 129 157 134 149 126 ...
##  $ hemoglobin         : num  12.9 12.7 15.8 14.7 12.5 16.2 17 15 13.7 16 ...
##  $ serum.creatinine   : num  0.7 0.6 1 1 0.6 1.2 0.7 1.3 0.8 0.8 ...
##  $ AST                : num  18 22 21 19 16 18 21 38 31 26 ...
##  $ ALT                : num  19 19 16 26 14 27 27 71 31 24 ...
##  $ Gtp                : num  27 18 22 18 22 33 39 111 14 63 ...
##  $ dental.caries      : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 2 1 1 1 ...
##  $ tartar             : Factor w/ 2 levels "N","Y": 2 2 1 2 1 2 2 2 1 1 ...
##  $ smoking            : Factor w/ 2 levels "0","1": 1 1 2 1 1 1 2 1 1 1 ...

3 Splitting dei dati

Innanzitutto si suddivida il dataset in training e test set. Il primo verrà usato per la fase di pre-processing e di scelta e valutazione dei modelli, il secondo per la classificazione delle nuove osservazioni e valutazione delle previsioni. Avendo una numerosità elevata di osservazioni si è scelto di effettuare una divisione 70% train - 30% test.

set.seed(123)
index <- sample(1:nrow(smoke),0.7*nrow(smoke)) 
train <- smoke[index,] #train: 38984 obs
test <- smoke[-index,] #test: 16708 obs

Dato che saranno valutate le performance di più modelli, si applica un’ulteriore divisione all’interno del set di training creando così un nuovo training set e un validation set.

sub.index <- sample(1:nrow(train),0.65*nrow(train))
sub.train <- train[sub.index,] #sub train: 25339 obs
validation <- train[-sub.index,] #validation: 13645 obs

Si verifica infine che le proporzioni delle due classi della variabile target siano le stesse su tutti i set di dati creati e che corrispondano a quelle del dataset di partenza (63% non fumatore - 37% fumatore).

round(prop.table(table(sub.train$smoking)),2)

## 
##    0    1 
## 0.63 0.37

round(prop.table(table(validation$smoking)),2)

## 
##    0    1 
## 0.63 0.37

round(prop.table(table(test$smoking)),2)

## 
##    0    1 
## 0.63 0.37

Tutto è correttamente proporzionato, con un leggero sbilanciamento delle classi. Si procede quindi con la fase di pre-processing.

4 Pre-processing

4.1 Analisi missing values e outliers

La seguente analisi esplorativa prenderà in considerazione solo il sub.train. Per prima cosa si ispeziona la presenza di valori mancanti.

sum(is.na(sub.train)) #sono presenti missing values?

## [1] 0

md.pattern(sub.train, rotate.names=TRUE)

##  /\     /\
## {  `---'  }
## {  O   O  }
## ==>  V <==  No need for mice. This data set is completely observed.
##  \  \|/  /
##   `-----'

##       gender age height.cm. weight.kg. waist.cm. hearing.left. hearing.right.
## 25339      1   1          1          1         1             1              1
##            0   0          0          0         0             0              0
##       systolic relaxation fasting.blood.sugar Cholesterol triglyceride HDL LDL
## 25339        1          1                   1           1            1   1   1
##              0          0                   0           0            0   0   0
##       hemoglobin serum.creatinine AST ALT Gtp dental.caries tartar smoking  
## 25339          1                1   1   1   1             1      1       1 0
##                0                0   0   0   0             0      0       0 0

sub.train.numeric <- sub.train%>% dplyr::select(-c(gender, hearing.left., hearing.right., dental.caries, tartar, smoking))
summary(sub.train.numeric) #statistiche principali per le variabili quantitative

##       age         height.cm.      weight.kg.       waist.cm.     
##  Min.   :20.0   Min.   :135.0   Min.   : 30.00   Min.   : 51.00  
##  1st Qu.:40.0   1st Qu.:160.0   1st Qu.: 55.00   1st Qu.: 76.00  
##  Median :40.0   Median :165.0   Median : 65.00   Median : 82.00  
##  Mean   :44.2   Mean   :164.7   Mean   : 65.86   Mean   : 82.03  
##  3rd Qu.:55.0   3rd Qu.:170.0   3rd Qu.: 75.00   3rd Qu.: 88.00  
##  Max.   :85.0   Max.   :190.0   Max.   :130.00   Max.   :128.00  
##     systolic       relaxation     fasting.blood.sugar  Cholesterol   
##  Min.   : 80.0   Min.   : 40.00   Min.   : 46.00      Min.   : 55.0  
##  1st Qu.:112.0   1st Qu.: 70.00   1st Qu.: 89.00      1st Qu.:172.0  
##  Median :120.0   Median : 76.00   Median : 96.00      Median :195.0  
##  Mean   :121.5   Mean   : 75.98   Mean   : 99.29      Mean   :196.8  
##  3rd Qu.:130.0   3rd Qu.: 82.00   3rd Qu.:104.00      3rd Qu.:219.0  
##  Max.   :213.0   Max.   :146.00   Max.   :505.00      Max.   :441.0  
##   triglyceride        HDL              LDL         hemoglobin   
##  Min.   : 11.0   Min.   :  4.00   Min.   :   1   Min.   : 4.90  
##  1st Qu.: 74.0   1st Qu.: 47.00   1st Qu.:  92   1st Qu.:13.60  
##  Median :107.0   Median : 55.00   Median : 113   Median :14.80  
##  Mean   :126.3   Mean   : 57.22   Mean   : 115   Mean   :14.62  
##  3rd Qu.:159.0   3rd Qu.: 65.00   3rd Qu.: 136   3rd Qu.:15.70  
##  Max.   :548.0   Max.   :618.00   Max.   :1660   Max.   :20.90  
##  serum.creatinine       AST               ALT               Gtp        
##  Min.   : 0.1000   Min.   :   6.00   Min.   :   2.00   Min.   :  3.00  
##  1st Qu.: 0.8000   1st Qu.:  19.00   1st Qu.:  15.00   1st Qu.: 17.00  
##  Median : 0.9000   Median :  23.00   Median :  21.00   Median : 25.00  
##  Mean   : 0.8839   Mean   :  26.28   Mean   :  27.17   Mean   : 39.93  
##  3rd Qu.: 1.0000   3rd Qu.:  28.00   3rd Qu.:  31.00   3rd Qu.: 43.00  
##  Max.   :10.0000   Max.   :1311.00   Max.   :2062.00   Max.   :999.00

Il dataset non contiene missing. E’ necessario controllare anche la presenza di eventuali outliers. Per fare ciò si riporta una tabella contenente massimo e minimo delle variabili quantitative e il range di variazione ottimale di esse (valori per una persona sana).

Train:

Variabile	Minimo	Massimo	Range di variazione ottimale
Systolic	80	213	110-140 mm/hg
Relaxation	40	146	70-90 mm/hg
Fasting blood sugar	46	505	70-140 mg/dl
Cholesterol	55	441	<200 mg/dl
Triglyceride	11	548	<200 mg/dl
HDL	4	618	>40 mg/dl
LDL	1	1660	<130 mg/dl
Hemoglobin	4.9	20.9	12-18 mg/dl
Serum creatinine	0.1	10	0.8-1.2 mg/dl
AST	6	1311	8-60 U/l
ALT	1	2062	7-55 U/l
γ-GTP	3	999	2-50 U/l

Essendo misurazioni mediche è normale che molte osservazioni siano fuori dal range di valori ottimale: ad esempio, un valore di 455 per il colesterolo è molto elevato e pericoloso per la salute, ma non per questo è da considerarsi un valore non plausibile. La presenza di outlier è quindi intrinseca nel tipo di variabili rilevate, inoltre la numerosità dei valori più estremi rispetto al totale delle osservazioni è irrisoria: queste considerazioni hanno portato alla conclusione di non rimuoverli, non causando particolari distorsioni sulle analisi.

4.2 Analisi delle correlazioni

correlazione <- cor(sub.train.numeric)
heatmaply::heatmaply_cor(correlazione, 
                         cellnote = round(correlazione,2), cellnote_textposition = "middle center",
                         cellnote_size = 8,
                         show_dendrogram = c(FALSE, FALSE))

Si notano alcune correlazioni significative:

• ALT - AST (83%)
• weight - waist (82%)
• relaxation - sistolic (76%)
• LDL - cholesterol (75%)

Le prime due coppie hanno una correlazione molto elevata, il che giustifica un’ulteriore selezione di variabili. In particolare, si mantengono AST (enzima più generico) e weight (waist fornisce un’informazione ridondante).

train <- train%>% dplyr::select(-c(waist.cm., ALT))
sub.train <- sub.train%>% dplyr::select(-c(waist.cm., ALT))
sub.train.numeric <- sub.train.numeric%>% dplyr::select(-c(waist.cm., ALT))
validation <- validation%>% dplyr::select(-c(waist.cm., ALT))
test <- test%>% dplyr::select(-c(waist.cm., ALT))

4.3 Distribuzione delle variabili

Si procede ora con una rappresentazione grafica delle esplicative attraverso dei boxplot condizionati alle classi della variabile target.

boxplot <-list()
variabili<- colnames(sub.train.numeric)
for (i in variabili){
  boxplot[[i]] <- ggplot(sub.train, aes_string(x = "smoking", y = i, col = "smoking", fill = "smoking")) + 
    geom_boxplot(alpha = 0.2) + 
    theme(legend.position = "none") + 
    scale_color_manual(values = c("blue", "red")) 
  scale_fill_manual(values = c("blue", "red"))
}
library(gridExtra)
grid.arrange(grobs=boxplot, nrow=4, ncol=4)

In alcuni boxplot si notano alcuni outlier: in HDL, ad esempio, c’è solo un valore completamente fuori dal range. Si tratta comunque di un osservazione su un sub.train con una numerosità di 25000, per questo non influenzerà le analisi e non verrà quindi eliminato. Si nota inoltre come alcune variabili abbiano una distribuzione fortemente asimmetrica, evidenziata anche dagli istogrammi di esse: per questo è opportuno effettuare una trasformazione logaritmica.

par(mfrow=c(1,2))

hist(sub.train$fasting.blood.sugar, main = 'fasting.blood.sugar')
sub.train$fasting.blood.sugar<-log(sub.train$fasting.blood.sugar)
hist(sub.train$fasting.blood.sugar, main = 'log(fasting.blood.sugar)')

hist(sub.train$triglyceride, main = 'triglyceride')
sub.train$triglyceride<-log(sub.train$triglyceride)
hist(sub.train$triglyceride, main = 'log(triglyceride)')

hist(sub.train$HDL, main = 'HDL')
sub.train$HDL<-log(sub.train$HDL)
hist(sub.train$HDL, main = 'log(HDL)')

hist(sub.train$LDL, main = 'LDL')
sub.train$LDL<-log(sub.train$LDL)
hist(sub.train$LDL, main = 'log(LDL)')

hist(sub.train$serum.creatinine, main = 'serum.creatinine')
sub.train$serum.creatinine<-log(sub.train$serum.creatinine)
hist(sub.train$serum.creatinine, main = 'log(serum.creatinine)')

hist(sub.train$AST, main = 'AST')
sub.train$AST<-log(sub.train$AST)
hist(sub.train$AST, main = 'log(AST)')

hist(sub.train$Gtp, main = 'Gtp')
sub.train$Gtp<-log(sub.train$Gtp)
hist(sub.train$Gtp, main = 'log(Gtp)')

par(mfrow=c(1,1))

Confrontando gli istogrammi delle variabili originali e delle trasformate logaritmiche si nota come la distribuzione sia ora più simmetrica.

La normalizzazione è altamente influenzata da quelli che potrebbero essere dei valori anomali nel sub.train, andando a genererare distorsioni nel validation e successivamente nel test, per questo motivo si è scelto di non utilizzarla. La scelta è ricaduta tra due tecniche, la prima è una standardizzazione e la seconda una tecnica costruita ad hoc andando a sottrarre al valore il primo quartile e dividendo il tutto per lo scarto interquartile. Si è preferito usare la prima: nonostante la media sia influenzata dagli outlier, la numerosità elevata mitiga questo effetto e si ottiene una variabile con distribuzione nota.

#Creazione matrice per tenerne traccia di min e max di ogni variabile 
matrix_indicators <- matrix(0, nrow=(dim(sub.train.numeric)[2]), ncol=2) 
colnames(matrix_indicators) <- c("media", "dv")

sub.train.numeric$fasting.blood.sugar<-(sub.train$fasting.blood.sugar)
sub.train.numeric$triglyceride<-(sub.train$triglyceride)
sub.train.numeric$HDL<-(sub.train$HDL)
sub.train.numeric$LDL<-(sub.train$LDL)
sub.train.numeric$serum.creatinine<-(sub.train$serum.creatinine)
sub.train.numeric$AST<-(sub.train$AST)
sub.train.numeric$Gtp<-(sub.train$Gtp)

#Calcolo media e deviazione standard 
for (i in 1:nrow(matrix_indicators)){
  matrix_indicators[i, "media"] <- mean(sub.train.numeric[,i])
  matrix_indicators[i, "dv"] <- sd(sub.train.numeric[,i])
}

#Standardizzazione
for (i in 1:nrow(matrix_indicators)){
  sub.train.numeric[, i] <- (sub.train.numeric[, i] - matrix_indicators[i, "media"])/(matrix_indicators[i, "dv"]) }

#Aggiornamento sub.train completo (con anche le qualitative)
sub.train <- cbind(sub.train.numeric, sub.train$gender, sub.train$hearing.left., sub.train$hearing.right.,
               sub.train$dental.caries, sub.train$tartar, sub.train$smoking)
colnames(sub.train)[15:20] <- c('gender','hearing.left.','hearing.right.','dental.caries','tartar','smoking')

4.4 Pca

Viene svolta la analisi per componenti principali al fine sia di individuare le variabili che più influenzano i pesi delle componenti, sia per graficare la distribuzione della classificazione con le due variabili più legate alle prime due componenti principali. Queste sono trigliceridi e Gtp: verranno tenute in considerazione solo a fini grafici, mentre i modelli verranno stimati con il set di variabili completo.

pca <- princomp(sub.train.numeric[,-c(1:3)])
summary(pca)

## Importance of components:
##                           Comp.1    Comp.2    Comp.3    Comp.4     Comp.5
## Standard deviation     1.7161246 1.3626498 1.1958172 1.0531326 1.03452262
## Proportion of Variance 0.2677454 0.1688080 0.1300032 0.1008302 0.09729812
## Cumulative Proportion  0.2677454 0.4365534 0.5665566 0.6673868 0.76468490
##                            Comp.6     Comp.7     Comp.8     Comp.9   Comp.10
## Standard deviation     0.90435635 0.76818257 0.72088586 0.61261831 0.4849610
## Proportion of Variance 0.07435388 0.05364798 0.04724518 0.03411964 0.0213815
## Cumulative Proportion  0.83903878 0.89268676 0.93993194 0.97405157 0.9954331
##                            Comp.11
## Standard deviation     0.224129940
## Proportion of Variance 0.004566928
## Cumulative Proportion  1.000000000

#le prime due componenti 
sort(pca$loadings[,1]) #Gtp

##                 HDL                 LDL         Cholesterol fasting.blood.sugar 
##         -0.24483128          0.09973841          0.15332011          0.23547218 
##    serum.creatinine                 AST            systolic          relaxation 
##          0.24401486          0.28763844          0.34607437          0.35833805 
##          hemoglobin        triglyceride                 Gtp 
##          0.35970767          0.38665206          0.42639597

sort(pca$loadings[,2]) #Cholesterol 

##         Cholesterol                 LDL                 HDL        triglyceride 
##         -0.69058207         -0.67570481         -0.18595675         -0.01900431 
##                 AST                 Gtp          relaxation          hemoglobin 
##          0.02604915          0.04414287          0.04468465          0.06289490 
##            systolic    serum.creatinine fasting.blood.sugar 
##          0.07491514          0.08022087          0.10464722

ggplot(sub.train, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point(alpha=0.2)+
  stat_ellipse(type="norm")+
  scale_colour_manual(values=c("darkgreen", "red"))+
  theme_bw()+
  labs(title="Classificazione originale")+
  xlab("trasformazione Gtp")+
  ylab("trasformazione Colesterolo")

Ogni modello verrà costruito sul sub train e testato nel validation, verrà creata una matrice che salva le seguenti metriche: F1 score, accuracy e recall(sensitivity). L’F1 score risulta particolarmente vantaggiosa in questo caso per due motivi. In primis è utile quando ci sono costi diversi di falsi positivi e falsi negativi: in questo caso, per un’assicurazione è peggio predire come non fumatori persone che in realtà lo sono, in quanto i costi calcolati sarebbero minore di quelli potenzialmente necessari. Inoltre le classi non sono perfettamente bilanciate (63%-37%), quindi limitarsi all’accuracy potrebbe dare risultati imprecisi.

matrice_ris<-matrix(ncol=4)
colnames(matrice_ris)<-c("modello", "F1", "accuracy", "recall")

L’intera analisi dei dati, data anche la numerosità del dataset, verrà effettuata valutando le prestazioni dei modelli costruiti sul sub.train nel validation. Si confronteranno successivamente i risultati ottenuti nel validation al fine di scegliere un unico modello da testare sul test set. L’obiettivo è simulare una situazione reale in cui il test set rappresenta i soggetti su cui usare il modello per predire la loro condizione di fumatori.

5 Modelli: training e validation

5.1 Alberi decisionali

Il primo modello testato è un albero di classificazione. I limiti di questo modello sono in primis a livello computazionale. Difatti vista la grande dimensionalità del dataset il numero di nodi che il modello va a creare è estremamente esiguo (3). Si è provato a rieseguire l’analisi sia riducendo la dimensionalità che usando le variabili senza le trasformazioni ma non sono stati ottenuti risultati migliori. Inoltre non è stata necessaria l’applicazione della potatura dell’albero in quanto già di dimensioni ridotte.

albero <- tree(smoking~., data=sub.train)
summary(albero)

## 
## Classification tree:
## tree(formula = smoking ~ ., data = sub.train)
## Variables actually used in tree construction:
## [1] "gender" "Gtp"   
## Number of terminal nodes:  3 
## Residual mean deviance:  0.9717 = 24620 / 25340 
## Misclassification error rate: 0.2891 = 7325 / 25339

plot(albero)
text(albero, pretty = 0)

L’albero ottenuto classifica inizialmente sul gender: se il genere è femminile la classificazione prevista è non fumatrici, se il genere è maschile invece viene effettuata una discriminazione in base al gtp. Come è immaginabile un modello così semplice risulta poco informativo, tanto che la misclassificazione è quasi del 29%.

albero

## node), split, n, deviance, yval, (yprob)
##       * denotes terminal node
## 
## 1) root 25339 33290 0 ( 0.63396 0.36604 )  
##   2) gender: F 9259  3063 0 ( 0.96079 0.03921 ) *
##   3) gender: M 16080 22100 1 ( 0.44577 0.55423 )  
##     6) Gtp < 0.606893 10756 14910 0 ( 0.50958 0.49042 ) *
##     7) Gtp > 0.606893 5324  6650 1 ( 0.31687 0.68313 ) *

Nello split basato sul Gtp la proporzione delle classi è quasi identica (55%-45%), questo rende la discriminazione legata al training set. Inoltre è interessante notare la proporzione di donne fumatrici (meno del 4%).

Prima di testare il modello sul validation, è opportuno e necessario effettuare su esso le trasformazioni svolte nel sub train.

summary(validation[,-20]) #no missing values

##  gender        age          height.cm.      weight.kg.     hearing.left.
##  F:4960   Min.   :20.00   Min.   :135.0   Min.   : 30.00   1:13292      
##  M:8685   1st Qu.:40.00   1st Qu.:160.0   1st Qu.: 55.00   2:  353      
##           Median :40.00   Median :165.0   Median : 65.00                
##           Mean   :44.14   Mean   :164.6   Mean   : 65.82                
##           3rd Qu.:55.00   3rd Qu.:170.0   3rd Qu.: 75.00                
##           Max.   :85.00   Max.   :190.0   Max.   :135.00                
##  hearing.right.    systolic       relaxation     fasting.blood.sugar
##  1:13286        Min.   : 79.0   Min.   : 46.00   Min.   : 46.0      
##  2:  359        1st Qu.:112.0   1st Qu.: 70.00   1st Qu.: 89.0      
##                 Median :120.0   Median : 76.00   Median : 96.0      
##                 Mean   :121.6   Mean   : 76.12   Mean   : 99.5      
##                 3rd Qu.:130.0   3rd Qu.: 82.00   3rd Qu.:104.0      
##                 Max.   :240.0   Max.   :146.00   Max.   :398.0      
##   Cholesterol     triglyceride        HDL              LDL        
##  Min.   : 77.0   Min.   : 16.0   Min.   : 18.00   Min.   :   1.0  
##  1st Qu.:172.0   1st Qu.: 74.0   1st Qu.: 47.00   1st Qu.:  91.0  
##  Median :195.0   Median :108.0   Median : 56.00   Median : 112.0  
##  Mean   :196.7   Mean   :126.6   Mean   : 57.41   Mean   : 114.9  
##  3rd Qu.:220.0   3rd Qu.:161.0   3rd Qu.: 66.00   3rd Qu.: 136.0  
##  Max.   :445.0   Max.   :999.0   Max.   :155.00   Max.   :1860.0  
##    hemoglobin    serum.creatinine       AST              Gtp        
##  Min.   : 4.90   Min.   : 0.1000   Min.   :  6.00   Min.   :  1.00  
##  1st Qu.:13.60   1st Qu.: 0.8000   1st Qu.: 19.00   1st Qu.: 17.00  
##  Median :14.80   Median : 0.9000   Median : 23.00   Median : 26.00  
##  Mean   :14.62   Mean   : 0.8854   Mean   : 26.23   Mean   : 40.21  
##  3rd Qu.:15.70   3rd Qu.: 1.0000   3rd Qu.: 29.00   3rd Qu.: 43.00  
##  Max.   :21.10   Max.   :11.6000   Max.   :981.00   Max.   :976.00  
##  dental.caries tartar  
##  0:10736       N:5986  
##  1: 2909       Y:7659  
##                        
##                        
##                        
## 

#Applicazione trasformate logaritmiche 
validation$fasting.blood.sugar<-log(validation$fasting.blood.sugar)
validation$triglyceride<-log(validation$triglyceride)
validation$HDL<-log(validation$HDL)
validation$LDL<-log(validation$LDL)
validation$serum.creatinine<-log(validation$serum.creatinine)
validation$AST<-log(validation$AST)
validation$Gtp<-log(validation$Gtp)

#Standardizzazione
validation.numeric <- validation[,-c(1,5,6,18,19,20)]
for (i in 1:nrow(matrix_indicators)){
  validation.numeric[, i] <- (validation.numeric[, i] - matrix_indicators[i, "media"])/(matrix_indicators[i, "dv"])
}
validation <- cbind(validation.numeric, validation$gender, validation$hearing.left., validation$hearing.right.,
               validation$dental.caries, validation$tartar, validation$smoking)
colnames(validation)[15:20] <- c('gender','hearing.left.','hearing.right.','dental.caries','tartar','smoking')

Si valuta l’efficienza sul validation.

tree.pred <- predict(albero , validation[,-20], type = "class")
cfalbero<-confusionMatrix(tree.pred , validation$smoking, positive = "1", mode="everything")

## Registered S3 methods overwritten by 'proxy':
##   method               from    
##   print.registry_field registry
##   print.registry_entry registry

cfalbero

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 7731 3050
##          1  909 1955
##                                           
##                Accuracy : 0.7099          
##                  95% CI : (0.7022, 0.7175)
##     No Information Rate : 0.6332          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.3136          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.3906          
##             Specificity : 0.8948          
##          Pos Pred Value : 0.6826          
##          Neg Pred Value : 0.7171          
##               Precision : 0.6826          
##                  Recall : 0.3906          
##                      F1 : 0.4969          
##              Prevalence : 0.3668          
##          Detection Rate : 0.1433          
##    Detection Prevalence : 0.2099          
##       Balanced Accuracy : 0.6427          
##                                           
##        'Positive' Class : 1               
## 

L’accuracy del 71% è fuorviante, infatti la classificazione dei non fumatori è efficiente (specificità quasi del 90%), ma facendo riferimento ai fumatori il modello risulta altamente inefficiente (sensibilità al di sotto del 40%).

#Salva i risultati ottenuti nella matrice
matrice_ris[1,]<-c("albero decionale", cfalbero$byClass["F1"], cfalbero$overall["Accuracy"],cfalbero$byClass["Sensitivity"])

Ora si provano a graficare i risultati ottenuti con la classificazione e le vere label, evidenziando in blu i punti misclassificati.

library(mclust)

## Warning: package 'mclust' was built under R version 4.2.3

## Package 'mclust' version 6.0.1
## Type 'citation("mclust")' for citing this R package in publications.

## 
## Attaching package: 'mclust'

## The following object is masked from 'package:purrr':
## 
##     map

misclassificatealbero<-(classError(tree.pred, validation$smoking)$misclassified)

plotv<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point()+
  scale_colour_manual(values=c("darkgreen", "red"))+
  theme_bw()+
  labs("Classificazione originale vere label")+
  theme_bw()
plotc<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point()+
  scale_colour_manual(values=c("darkgreen", "red"), name="classificazione")+
  geom_point(data=validation[misclassificatealbero,],  col="blue", alpha=0.5)+
  theme_bw()+
  labs("Classificazione secondo il metodo tree")+
  theme_bw()
grid.arrange(plotv, plotc, ncol = 2,
             top = "Etichette vere vs classificazione con un albero")

Il blu indica una classificazione errata. Nel grafico si nota come l’albero abbia una grossa difficoltà a classificare principalmente le unità “centrali”.

5.2 LDA e QDA

Il modello più restrittivo in termini di ipotesi sulle covariate è quello basato su LDA/QDA. Questi metodi impongono in primis che sia verificata la distribuzione normale delle variabili. Per avere una visione immediata vengono qui presentati i QQplot condizionati alla classe; ad essi viene sovraimposta la linea dei quantili ‘teorici’ in modo da avere un’indicazione dello scostamento di quelli della variabile in esame.

Da questa analisi sono escluse le variabili dicotomiche, age, height e weight. Osservando i valori di queste ultime 3 si nota infatti come assumano valori ad intervalli di 5 unità: le rilevazioni sono quindi state arrotondate. Si potrebbe infatti pensare di ricodificarle utilizzando degli intervalli e renderle variabili factor. Al di fuori di esse rimane comunque un ampio set di variabili, quindi si procede con esse.

#QQplot condizionati alla classe smoking = 1 (fumatori)
par(mfrow = c(2,3))
for(i in variabili[4:9]) {
  qqnorm(sub.train[sub.train$smoking == 1, i], main = i)
  qqline(sub.train[sub.train$smoking == 1, i], col = 2)
}

par(mfrow = c(2,3))
for(i in variabili[10:14]) {
  qqnorm(sub.train[sub.train$smoking == 1, i], main = i)
  qqline(sub.train[sub.train$smoking == 1, i], col = 2)
}

#QQplot condizionati alla classe smoking = 0 (non fumatori)
par(mfrow = c(2,3))

for(i in variabili[4:9]) {
  qqnorm(sub.train[sub.train$smoking == 0, i], main = i)
  qqline(sub.train[sub.train$smoking == 0, i], col = 2)
}

par(mfrow = c(2,3))
for(i in variabili[10:14]) {
  qqnorm(sub.train[sub.train$smoking == 0, i], main = i)
  qqline(sub.train[sub.train$smoking == 0, i], col = 2)
}

Già dai QQplot si nota come nessuna delle variabili abbia un andamento normale per entrambe le classi. Per avere un’ulteriore conferma si costruiscano le curve di densità condizionate.

plot_density <- list()
for(i in variabili[4:14]){
  plot_density[[i]] <- ggplot(sub.train, aes_string(x = i, y = "..density..", col = "smoking")) + 
    geom_density(aes(y = ..density..)) + 
    scale_color_manual(values = c("blue", "red")) + 
    theme(legend.position = "none")
}
do.call(grid.arrange, c(plot_density, nrow = 4))

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

I problemi sono soprattutto sulle code: questo è causato dai numerosi valori estremi che caratterizzano le variabili. Un’ulteriore conferma sarebbe fornita dal test per la normalità di Shapiro-Wilk, non applicabile in quest’analisi a causa della numerosità delle osservazioni molto maggiore rispetto a quella richiesta per lo svolgimento del test (dal punto di vista computazionale).

Dall’analisi fatta a livello qualitativo sui grafici delle densità, alcune variabili come colesterolo o HDL hanno una distribuzione normale. Il problema risiede nella perfetta sovrapposizione delle due curve di densità condizionate, che impediscono ai modelli di identificare in maniera efficace le due classi. E’ stata comunque svolta un’analisi per un’ulteriore conferma, la quale ha dato pessimi risultati. Le uniche variabili che sembrano essere distinte sono i trigliceridi, emoglobina e Gtp ma non presentano una distribuzione normale, potrebbero essere misture di misture. Si sceglie di procedere con metodi ritenuti più opportuni.

5.3 Regressione logistica

Si applica ora la regressione logistica. Prima si stima un modello con tutte le esplicative e successivamente, per cercare di ridurre la dimensionalità, si usa una selezione stepwise basata sul criterio AIC (the lower the better) in modo da rimuovere le variabili non significative.

logit <- glm(smoking ~., data = sub.train, family = binomial)
summary(logit) #AIC: 23405

## 
## Call:
## glm(formula = smoking ~ ., family = binomial, data = sub.train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5620  -0.8798  -0.2321   0.8924   3.0405  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -3.213295   0.067805 -47.391  < 2e-16 ***
## age                 -0.025855   0.018746  -1.379  0.16782    
## height.cm.           0.211427   0.029391   7.194 6.31e-13 ***
## weight.kg.          -0.220915   0.024554  -8.997  < 2e-16 ***
## systolic            -0.197098   0.025913  -7.606 2.83e-14 ***
## relaxation           0.045124   0.025163   1.793  0.07293 .  
## fasting.blood.sugar  0.026668   0.017003   1.568  0.11678    
## Cholesterol         -0.001838   0.050627  -0.036  0.97104    
## triglyceride         0.276122   0.028689   9.625  < 2e-16 ***
## HDL                 -0.026433   0.027175  -0.973  0.33071    
## LDL                 -0.119755   0.045465  -2.634  0.00844 ** 
## hemoglobin           0.171546   0.025415   6.750 1.48e-11 ***
## serum.creatinine    -0.200441   0.021282  -9.418  < 2e-16 ***
## AST                 -0.226599   0.018803 -12.051  < 2e-16 ***
## Gtp                  0.537373   0.022545  23.835  < 2e-16 ***
## genderM              2.992413   0.079223  37.772  < 2e-16 ***
## hearing.left.2      -0.115485   0.120547  -0.958  0.33806    
## hearing.right.2      0.092579   0.119591   0.774  0.43886    
## dental.caries1       0.310428   0.038648   8.032 9.57e-16 ***
## tartarY              0.306661   0.032797   9.350  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 33286  on 25338  degrees of freedom
## Residual deviance: 23365  on 25319  degrees of freedom
## AIC: 23405
## 
## Number of Fisher Scoring iterations: 6

step.logit <- stepAIC(logit, direction = "both", trace = FALSE)
summary(step.logit) #AIC: 23400

## 
## Call:
## glm(formula = smoking ~ height.cm. + weight.kg. + systolic + 
##     relaxation + triglyceride + LDL + hemoglobin + serum.creatinine + 
##     AST + Gtp + gender + dental.caries + tartar, family = binomial, 
##     data = sub.train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5896  -0.8816  -0.2323   0.8921   3.0394  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -3.22105    0.06733 -47.838  < 2e-16 ***
## height.cm.        0.21405    0.02795   7.659 1.88e-14 ***
## weight.kg.       -0.20590    0.02353  -8.752  < 2e-16 ***
## systolic         -0.20045    0.02564  -7.818 5.36e-15 ***
## relaxation        0.04526    0.02512   1.802   0.0716 .  
## triglyceride      0.28822    0.01847  15.601  < 2e-16 ***
## LDL              -0.12326    0.01634  -7.543 4.60e-14 ***
## hemoglobin        0.17642    0.02494   7.075 1.50e-12 ***
## serum.creatinine -0.20238    0.02122  -9.537  < 2e-16 ***
## AST              -0.22901    0.01871 -12.238  < 2e-16 ***
## Gtp               0.53476    0.02208  24.218  < 2e-16 ***
## genderM           3.00038    0.07799  38.473  < 2e-16 ***
## dental.caries1    0.31451    0.03846   8.179 2.87e-16 ***
## tartarY           0.30991    0.03275   9.463  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 33286  on 25338  degrees of freedom
## Residual deviance: 23372  on 25325  degrees of freedom
## AIC: 23400
## 
## Number of Fisher Scoring iterations: 6

Le variabili significative mantenute dalla selezione stepwise sono comunque molto numerose. Si procede analizzando la presenza di punti influenti, eliminandoli e ripetendo la stima del modello per vedere se ci sono dei miglioramenti.

influencePlot(step.logit)

##         StudRes          Hat        CookD
## 20340 -1.675502 8.783723e-03 0.0019259620
## 4225  -2.468969 2.086182e-03 0.0029393022
## 11570  1.682581 8.739313e-03 0.0019461947
## 1577   3.040595 6.957062e-05 0.0004990126
## 17822 -2.598540 1.683238e-03 0.0033280060
## 16607  3.004504 8.576728e-05 0.0005507912

sub.train2 <- sub.train[-c(20340,4225,11570,1577,17822,16607),]
logit2 <- glm(smoking ~., data = sub.train2, family = binomial)
summary(logit2) #AIC: 23402

## 
## Call:
## glm(formula = smoking ~ ., family = binomial, data = sub.train2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5612  -0.8802  -0.2320   0.8925   3.0398  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         -3.212615   0.067802 -47.382  < 2e-16 ***
## age                 -0.025509   0.018747  -1.361  0.17362    
## height.cm.           0.211562   0.029392   7.198 6.11e-13 ***
## weight.kg.          -0.220420   0.024556  -8.976  < 2e-16 ***
## systolic            -0.197443   0.025917  -7.618 2.57e-14 ***
## relaxation           0.045226   0.025164   1.797  0.07230 .  
## fasting.blood.sugar  0.026474   0.017007   1.557  0.11955    
## Cholesterol          0.001999   0.050975   0.039  0.96871    
## triglyceride         0.274407   0.028788   9.532  < 2e-16 ***
## HDL                 -0.028021   0.027269  -1.028  0.30415    
## LDL                 -0.122991   0.045757  -2.688  0.00719 ** 
## hemoglobin           0.171089   0.025417   6.731 1.68e-11 ***
## serum.creatinine    -0.200285   0.021282  -9.411  < 2e-16 ***
## AST                 -0.226378   0.018804 -12.039  < 2e-16 ***
## Gtp                  0.537272   0.022544  23.832  < 2e-16 ***
## genderM              2.992071   0.079222  37.768  < 2e-16 ***
## hearing.left.2      -0.115652   0.120543  -0.959  0.33735    
## hearing.right.2      0.092334   0.119589   0.772  0.44005    
## dental.caries1       0.310326   0.038650   8.029 9.82e-16 ***
## tartarY              0.306032   0.032799   9.331  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 33277  on 25332  degrees of freedom
## Residual deviance: 23362  on 25313  degrees of freedom
## AIC: 23402
## 
## Number of Fisher Scoring iterations: 6

step.logit2 <- stepAIC(logit2, direction = "both", trace = FALSE)
summary(step.logit2) #AIC: 23396

## 
## Call:
## glm(formula = smoking ~ height.cm. + weight.kg. + systolic + 
##     relaxation + triglyceride + LDL + hemoglobin + serum.creatinine + 
##     AST + Gtp + gender + dental.caries + tartar, family = binomial, 
##     data = sub.train2)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.5886  -0.8817  -0.2324   0.8921   3.0390  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)      -3.22038    0.06733 -47.831  < 2e-16 ***
## height.cm.        0.21398    0.02795   7.656 1.92e-14 ***
## weight.kg.       -0.20542    0.02353  -8.730  < 2e-16 ***
## systolic         -0.20081    0.02564  -7.831 4.83e-15 ***
## relaxation        0.04541    0.02512   1.808   0.0706 .  
## triglyceride      0.28811    0.01847  15.596  < 2e-16 ***
## LDL              -0.12326    0.01634  -7.543 4.60e-14 ***
## hemoglobin        0.17586    0.02494   7.052 1.77e-12 ***
## serum.creatinine -0.20214    0.02122  -9.525  < 2e-16 ***
## AST              -0.22868    0.01871 -12.220  < 2e-16 ***
## Gtp               0.53467    0.02208  24.215  < 2e-16 ***
## genderM           3.00001    0.07798  38.470  < 2e-16 ***
## dental.caries1    0.31436    0.03846   8.174 2.98e-16 ***
## tartarY           0.30925    0.03275   9.442  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 33277  on 25332  degrees of freedom
## Residual deviance: 23368  on 25319  degrees of freedom
## AIC: 23396
## 
## Number of Fisher Scoring iterations: 6

Il modello risultante comprende le seguenti variabili: height, weight, systolic, relaxation, triglyceride, LDL, hemoglobin, serum.creatinine, AST, Gtp, gender, dental.caries, tartar.

La regressione logistica assume che ci sia una relazione lineare tra le variabili esplicative e il logaritmo dell’odds ratio: è quindi bene testare ciò attraverso un’analisi grafica.

#Calcolo fitted values (focus solo sulle variabili quantitative)
probabilities <- predict(step.logit2, type = "response") #posterior probabilities
variabili.new <- variabili[-(6:7)]
#Costruzione degli odds ratio 
supp <- sub.train2[,variabili.new]
supp <- supp %>%
  mutate(logit = log(probabilities/(1-probabilities))) %>%
  gather(key = "variabili.new", value = "predictor.value", -logit)

#Costruzione dei grafici
ggplot(supp, aes(logit, predictor.value))+
  geom_point(size = 0.5, alpha = 0.5) +
  geom_smooth(method = "loess") + 
  theme_bw() + 
  facet_wrap(~variabili[-(6:7)], scales = "free_y")

## `geom_smooth()` using formula = 'y ~ x'

Come si nota la linearità è in parte rispettata ma l’alta variabilità dei punti non permette un’investigazione più precisa. Inoltre c’è un’ulteriore evidenza del fatto che age, height e weight andrebbero trasformate in dummy individuando più intervalli a cui assegnare le osservazioni.

predict.logit <- predict(step.logit2, validation[, -20], type = "response") #posterior probabilities
predict.logit.class <- ifelse(predict.logit > 0.5, 1, 0) #suddivide le post prob in due classi (0-1) con la soglia decisionale pari a 0.5 
conf <- confusionMatrix(factor(validation[, 20]), factor(as.vector(predict.logit.class)), positive = '1', mode='everything')
conf

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 6745 1895
##          1 1505 3500
##                                           
##                Accuracy : 0.7508          
##                  95% CI : (0.7435, 0.7581)
##     No Information Rate : 0.6046          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.4722          
##                                           
##  Mcnemar's Test P-Value : 2.536e-11       
##                                           
##             Sensitivity : 0.6487          
##             Specificity : 0.8176          
##          Pos Pred Value : 0.6993          
##          Neg Pred Value : 0.7807          
##               Precision : 0.6993          
##                  Recall : 0.6487          
##                      F1 : 0.6731          
##              Prevalence : 0.3954          
##          Detection Rate : 0.2565          
##    Detection Prevalence : 0.3668          
##       Balanced Accuracy : 0.7332          
##                                           
##        'Positive' Class : 1               
## 

La sensitività è meno del 65%, l’F1 score è circa il 67%. Perciò si provi a modificare la soglia per l’assegnazione al fine di migliorare la classificazione dei fumatori.

predict.logit.class2 <- ifelse(predict.logit > 0.3, 1, 0) 
conf1 <- confusionMatrix(factor(as.vector(predict.logit.class2)),factor(validation[, 20]),positive='1',
                         mode='everything')
conf1

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 5080  334
##          1 3560 4671
##                                          
##                Accuracy : 0.7146         
##                  95% CI : (0.707, 0.7222)
##     No Information Rate : 0.6332         
##     P-Value [Acc > NIR] : < 2.2e-16      
##                                          
##                   Kappa : 0.459          
##                                          
##  Mcnemar's Test P-Value : < 2.2e-16      
##                                          
##             Sensitivity : 0.9333         
##             Specificity : 0.5880         
##          Pos Pred Value : 0.5675         
##          Neg Pred Value : 0.9383         
##               Precision : 0.5675         
##                  Recall : 0.9333         
##                      F1 : 0.7058         
##              Prevalence : 0.3668         
##          Detection Rate : 0.3423         
##    Detection Prevalence : 0.6032         
##       Balanced Accuracy : 0.7606         
##                                          
##        'Positive' Class : 1              
## 

Il valore della recall è aumentato nettamente, arrivando al 93%. Anche l’F1 score è migliorato. Diminuendo ulteriormente la soglia migliora la recall ma peggiora la precision, quindi si sceglie di mantenere la soglia a 0.3.

matrice_ris<-rbind(matrice_ris,c("regressione logistica",conf1$byClass["F1"],conf1$overall["Accuracy"],conf1$byClass["Sensitivity"]))

Si graficano quindi i risultati ottenuti.

misclassificatereg<-(classError(predict.logit.class2, validation$smoking)$misclassified)

plotv<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point()+
  scale_colour_manual(values=c("darkgreen", "red"))+
  theme_bw()+
  labs("Classificazione originale vere label")+
  theme_bw()
plotc<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point()+
  scale_colour_manual(values=c("darkgreen", "red"), name="classificazione")+
  geom_point(data=validation[misclassificatealbero,],  col="blue", alpha=0.5)+
  theme_bw()+
  labs("Classificazione secondo la regressione logistica")+
  theme_bw()
grid.arrange(plotv, plotc, ncol = 2,
             top = "Etichette vere vs classificazione con regressione logistica")

Si valuti ora la curva ROC e si calcoli l’AUC.

roc.logit <- prediction(predict.logit, validation[, 20])
performance.logit <- performance(roc.logit,"tpr","fpr")
auc.logit <- performance(roc.logit, measure = "auc")@y.values
auc.logit

## [[1]]
## [1] 0.8351383

plot(performance.logit, colworize = TRUE, main = "Regressione Logistica")

Il modello di regressione logistica non si adatta bene ai dati per il problema già citato della linearità. Inoltre un problema riscontrato è dovuto alla dimensionalità. Sono infatti state effettuate più prove per ridurre il numero di variabili come: analisi delle componenti principali, eliminazione di variabili più specifiche (es. rimuovere HDL e LDL e mantenere Cholesterol). Questi tentativi hanno però condotto a modelli con un AIC maggiore e senza miglioramenti in termini di performance: questo ha portato alla decisione di quasi tutte le variabili, a discapito di un modello molto ricco.

5.4 KNN

Si sceglie ora di valutare un modello non parametrico basato sulle distanze: per questo motivo verrano usate solo le variabili quantitative già opportunamente trasformate. La funzione KNN chiede in input la risposta e le covariate in due dataset distinti, sia per il train che per il validation.

X_train_picc<-sub.train%>%dplyr::select(-c(smoking,gender,hearing.left., hearing.right., dental.caries, tartar))#SOLO QUANTITATIVE
Y_train_picc<-sub.train$smoking
#validation
X_validation<-validation%>%dplyr::select(-c(smoking,gender,hearing.left., hearing.right., dental.caries, tartar))
Y_validation<-validation$smoking

Si prova ora a trovare il numero di k vicini ottimale andando a confrontare i diversi valori dell’F1 per ciascun k.

k_to_try = 1:100
err_k = rep(x = 0, times = length(k_to_try))
for (i in seq_along(k_to_try)) {
  pred = knn(train = X_train_picc,
             test = X_validation, cl = Y_train_picc,
             k = k_to_try[i])
  err_k[i] = confusionMatrix(pred, Y_validation, positive='1')$byClass["F1"]}

L’esecuzione del ciclo risulta esssere onerosa computazionalmente ma i vantaggi operativi sono notevoli. Tramite un plot dei risultati ottenuti con l’esecuzione del ciclo for, si valutano graficamente le performance del modello al variare del parametro di tuning k.

F1_scores<-matrix(nrow=length(k_to_try), ncol=2)
F1_scores[,1]<-1:100
F1_scores[,2]<-err_k
colnames(F1_scores)<-c("k", "F1")
F1gg<-as.data.frame(F1_scores)
ggplot(data=F1gg, aes(x=k, y=F1))+
  geom_point(col="blue")+
  geom_line(lty="dotted")+
  theme_bw()+
  geom_text(data=F1gg[91,], aes(label=round(F1,2)), vjust=-0.8)+
  geom_hline(yintercept =F1gg[91,2], col="red", alpha=0.4)+
  scale_x_continuous(breaks = c(1,seq(5,100, by=5), 91) )+
  labs("f1 score in base al k")

Il miglior k è risultato essere quello uguale a 91. Si ristima il modello sul sub train valutandone la performance nel validation.

#miglior k
best_k<-knn(train = X_train_picc,
            test = X_validation, cl = Y_train_picc,
            k = 91, prob=T)

conf2 <- confusionMatrix(best_k, Y_validation, positive='1', mode='everything')
conf2

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 6864 1773
##          1 1776 3232
##                                           
##                Accuracy : 0.7399          
##                  95% CI : (0.7325, 0.7473)
##     No Information Rate : 0.6332          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.4401          
##                                           
##  Mcnemar's Test P-Value : 0.9732          
##                                           
##             Sensitivity : 0.6458          
##             Specificity : 0.7944          
##          Pos Pred Value : 0.6454          
##          Neg Pred Value : 0.7947          
##               Precision : 0.6454          
##                  Recall : 0.6458          
##                      F1 : 0.6456          
##              Prevalence : 0.3668          
##          Detection Rate : 0.2369          
##    Detection Prevalence : 0.3670          
##       Balanced Accuracy : 0.7201          
##                                           
##        'Positive' Class : 1               
## 

Il modello ha una sensitività e un F1 score non ottimali, rispettivamente intorno al 65%. Questo è dovuto alla struttura dei cluster in cui i punti sono molto sovrapposti tra le due classi. Una soluzione è quella della modifica delle soglie per far si che l’assegnazione avvenga in maniera più coerente rispetto alla reale proporzione delle classi.

prob.knn <- attributes(best_k)$prob
prop_1<-ifelse(best_k == "0", 1-prob.knn, prob.knn)
classi_soglia_knn<-ifelse(prop_1>0.4, 1, 0)
conf3<-confusionMatrix(as.factor(classi_soglia_knn), validation$smoking, positive = "1", mode="everything")
conf3

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 5676  804
##          1 2964 4201
##                                           
##                Accuracy : 0.7239          
##                  95% CI : (0.7163, 0.7313)
##     No Information Rate : 0.6332          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.455           
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.8394          
##             Specificity : 0.6569          
##          Pos Pred Value : 0.5863          
##          Neg Pred Value : 0.8759          
##               Precision : 0.5863          
##                  Recall : 0.8394          
##                      F1 : 0.6904          
##              Prevalence : 0.3668          
##          Detection Rate : 0.3079          
##    Detection Prevalence : 0.5251          
##       Balanced Accuracy : 0.7482          
##                                           
##        'Positive' Class : 1               
## 

La sensibilità, abbassando la soglia a 0.4, è aumentata fino a quasi l’84% e l’F1 score è diventato del 69%.

library(mclust)
misclassificateknn<-(classError(classi_soglia_knn, validation$smoking)$misclassified)

plotv<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point()+
  scale_colour_manual(values=c("darkgreen", "red"))+
  theme_bw()+
  labs("Classificazione originale vere label")+
  theme_bw()
plotc<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point()+
  scale_colour_manual(values=c("darkgreen", "red"), name="classificazione")+
  geom_point(data=validation[misclassificateknn,],  col="blue", alpha=0.5)+
  theme_bw()+
  labs("Classificazione secondo il metodo knn")+
  theme_bw()
grid.arrange(plotv, plotc, ncol = 2,
             top = "Etichette vere vs classificazione con 91-nn")

matrice_ris<-rbind(matrice_ris,c("91-NN",conf3$byClass["F1"],conf3$overall["Accuracy"],conf3$byClass["Sensitivity"]))

Si valuti ora la curva ROC e si calcoli l’AUC.

prob.knn <- attributes(best_k)$prob
prob.knn <- 2*ifelse(best_k == "0", 1-prob.knn, prob.knn) - 1
roc.knn <- prediction(prob.knn, validation[, 20])
performance.knn <- performance(roc.knn,"tpr","fpr")
auc.knn <- performance(roc.knn, measure = "auc")@y.values
auc.knn

## [[1]]
## [1] 0.8188348

plot(performance.knn, colworize = TRUE, main = "91-NN")

5.5 Random Forest

Il metodo Random Forest si basa sull’idea di combinare i risultati di più alberi decisionali.

Nella prima fase dell’algoritmo vengono addestrati più alberi decisionali a partire da campioni casuali estratti con reinserimento dal dataset di partenza (più training set). In questo modo si riduce l’overfitting: ciascun albero è legato ai suoi dati di training, ma la foresta non lo è.

Per quanto riguarda la costruzione di ciascun albero le features da valutare in ogni nodo sono selezionate casualmente: in questo modo gli alberi saranno meno correlati tra di loro e più diversi, riducendo la varianza e aumentando le prestazioni del modello.

Quando si considera una nuova osservazione basta ottenere la previsione della classe da ogni albero decisionale: la previsione finale è quella più frequente guardando all’intera foresta.

Dato che Random Forest non necessita di particolari ipotesi sulle variabili sono state fatte più prove sia usando il dataset originale che quello trasformato (trasformazioni logaritmiche e standardizzazione). I risultati sono pressochè identici, quindi per coerenza procediamo con le variabili standardizzate.

set.seed(123)
rf <- randomForest(smoking ~., data=sub.train, 
                   importance=TRUE, #l'argomento importance = TRUE fa sì che venga valutata l'importanza relativa di ciascun esplicativa
                   cutoff=c(0.7,0.3)) #regola la soglia decisionale

#valori previsti sul validation
rf.predict1 <- predict(rf, validation[,-20], type='prob')
rf.predict1.class <- predict(rf, validation[,-20])
conf4 <- confusionMatrix(rf.predict1.class, validation$smoking, positive='1', mode='everything')
conf4

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 5660  248
##          1 2980 4757
##                                           
##                Accuracy : 0.7634          
##                  95% CI : (0.7562, 0.7705)
##     No Information Rate : 0.6332          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.5432          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.9504          
##             Specificity : 0.6551          
##          Pos Pred Value : 0.6148          
##          Neg Pred Value : 0.9580          
##               Precision : 0.6148          
##                  Recall : 0.9504          
##                      F1 : 0.7467          
##              Prevalence : 0.3668          
##          Detection Rate : 0.3486          
##    Detection Prevalence : 0.5670          
##       Balanced Accuracy : 0.8028          
##                                           
##        'Positive' Class : 1               
## 

Il valore della sensitivity è ottimo: 95%. I risultati sono molto buoni anche in termini di F1 score. Il fatto di aver spostato la soglia è andato a discapito della specificity che però, per questa analisi, è meno rilevante.

misclassificaterf<-(classError(rf.predict1.class, validation$smoking)$misclassified)

plotv<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point()+
  scale_colour_manual(values=c("darkgreen", "red"))+
  theme_bw()+
  labs("Classificazione originale vere label")+
  theme_bw()
plotc<-ggplot(validation, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point()+
  scale_colour_manual(values=c("darkgreen", "red"), name="classificazione")+
  geom_point(data=validation[misclassificaterf,],  col="blue", alpha=0.5)+
  theme_bw()+
  labs("Classificazione secondo il metodo random forest")+
  theme_bw()
grid.arrange(plotv, plotc, ncol = 2,
             top = "Etichette vere vs classificazione con random forest")

matrice_ris<-rbind(matrice_ris,c("random forest",conf4$byClass["F1"],conf4$overall["Accuracy"],conf4$byClass["Sensitivity"]))

#curva ROC e AUC
roc.rf <- prediction(rf.predict1[,2], validation[, 20])
performance.rf <- performance(roc.rf,"tpr","fpr")
auc.rf <- performance(roc.rf, measure = "auc")@y.values
auc.rf

## [[1]]
## [1] 0.8882169

plot(performance.rf, colworize = TRUE, main = "Random Forest")

Il coefficiente AUC è di 0.89: si tratta di un ottimo risultato, considerando che un classificatore perfetto ha AUC=1.

Si proceda ora andando a investigare l’importanza relativa di ciascuna delle variabili esplicative:

importance <- importance(rf)
var.imp <- tibble(variabile=row.names(importance),
                  importanza = importance[,'MeanDecreaseGini'])
ggplot(var.imp, aes(x=reorder(variabile,desc(importanza)), y=importanza, fill=importanza))+
  geom_col()+
  labs(x="Predittori", y ="Importanza", title="Importanza dei predittori")+
  theme(axis.text.x = element_text(angle=90))

L’importanza delle covariate è tanto maggiore quanto è il valore dell’indice di impurità di Gini. Le variabili più importanti per la classificazione sono genere, gtp, emoglobina e in generale le variabili legate al colesterolo. Minor importanza hanno quelle legate all’udito e all’igiene dentale.

6 Valutazione sul test

Si confrontino innanzitutto i risultati ottenuti dai modelli precedenti in termini di metriche.

matrice_ris<-as.data.frame(matrice_ris)
matrice_ris$F1<-as.numeric(matrice_ris$F1)
matrice_ris$accuracy<-as.numeric(matrice_ris$accuracy)
matrice_ris$recall<-as.numeric(matrice_ris$recall)
matrice_ris

##                 modello        F1  accuracy    recall
## 1      albero decionale 0.4968865 0.7098571 0.3906094
## 2 regressione logistica 0.7058024 0.7146207 0.9332667
## 3                 91-NN 0.6903862 0.7238549 0.8393606
## 4         random forest 0.7466646 0.7634298 0.9504496

plot.a<-ggplot(matrice_ris, aes(x=modello,y=F1,col=modello))+
  geom_point()+
  geom_text(aes(label=round(F1,4)),col='black',vjust=1.5,size=3)+
  scale_y_continuous(limits=c(0,1))+
  theme_bw()+
  xlab('')+
  theme(legend.position = 'top')
plot.b<-ggplot(matrice_ris, aes(x=modello,y=accuracy,col=modello))+
  geom_point()+
  geom_text(aes(label=round(accuracy,4)),col='black',vjust=1.5,size=3)+
  scale_y_continuous(limits=c(0,1))+
  theme_bw()+
  xlab('')+
  theme(legend.position = 'none')
plot.c<-ggplot(matrice_ris, aes(x=modello,y=recall,col=modello))+
  geom_point()+
  geom_text(aes(label=round(recall,4)),col='black',vjust=1.5,size=3)+
  scale_y_continuous(limits=c(0,1))+
  theme_bw()+
  theme(legend.position = 'none')

grid.arrange(plot.a,plot.b,plot.c, nrow = 3, top='Confronto metriche sui vari modelli',heights=c(21.5,15,15))

Come si nota dal grafico il modello che performa meglio in termini di tutte le metriche è random forest. Il peggiore invece, come si era già riscontrato, è l’albero decisionale. L’accuracy non è la metrica più adeguata per questo problema, infatti dà uguale peso agli errori di classificazione. In questo caso invece è bene minimizzare il tasso di falsi negativi (non fumatori) e quindi massimizzare la recall: questa metrica è infatti molto più utile agli scopi di questa analisi ed è anche quella che discrimina maggiormente tra i metodi, avvalorando il fatto che random forest sia quello migliore. L’F1 combina recall e precisione, che misura la qualità della classificazione della classe di interesse (fumatori). Anche questa assume il valore più elevato con il metodo random forest.

Il fatto che quest’ultimo sia il classificatore migliore è visibile anche confrontando le curve di ROC.

par(mfrow=c(1,3))
plot(performance.logit, colworize = TRUE, main = "Regressione Logistica")
plot(performance.knn, colworize = TRUE, main = "91-KNN")
plot(performance.rf, colworize = TRUE, main = "Random Forest")

par(mfrow=c(1,1))

Si va quindi a valutare la performance di random forest sul test set. In primis si uniscano sub.train e validation e si trasformi il test set.

summary(test[,-20]) #no missing values

##  gender         age          height.cm.      weight.kg.     hearing.left.
##  F: 6072   Min.   :20.00   Min.   :130.0   Min.   : 35.00   1:16280      
##  M:10636   1st Qu.:40.00   1st Qu.:160.0   1st Qu.: 55.00   2:  428      
##            Median :40.00   Median :165.0   Median : 65.00                
##            Mean   :44.19   Mean   :164.6   Mean   : 65.92                
##            3rd Qu.:55.00   3rd Qu.:170.0   3rd Qu.: 75.00                
##            Max.   :85.00   Max.   :190.0   Max.   :125.00                
##  hearing.right.    systolic       relaxation     fasting.blood.sugar
##  1:16271        Min.   : 71.0   Min.   : 40.00   Min.   : 48.0      
##  2:  437        1st Qu.:112.0   1st Qu.: 70.00   1st Qu.: 89.0      
##                 Median :120.0   Median : 76.00   Median : 96.0      
##                 Mean   :121.4   Mean   : 75.96   Mean   : 99.2      
##                 3rd Qu.:130.0   3rd Qu.: 82.00   3rd Qu.:104.0      
##                 Max.   :220.0   Max.   :146.00   Max.   :475.0      
##   Cholesterol     triglyceride        HDL             LDL         hemoglobin   
##  Min.   : 77.0   Min.   :  8.0   Min.   :  4.0   Min.   :   1   Min.   : 5.00  
##  1st Qu.:173.0   1st Qu.: 75.0   1st Qu.: 47.0   1st Qu.:  92   1st Qu.:13.60  
##  Median :195.0   Median :109.0   Median : 55.0   Median : 113   Median :14.80  
##  Mean   :197.3   Mean   :127.2   Mean   : 57.3   Mean   : 115   Mean   :14.63  
##  3rd Qu.:220.0   3rd Qu.:160.0   3rd Qu.: 66.0   3rd Qu.: 136   3rd Qu.:15.80  
##  Max.   :393.0   Max.   :399.0   Max.   :359.0   Max.   :1660   Max.   :20.90  
##  serum.creatinine       AST           Gtp         dental.caries tartar  
##  Min.   : 0.1000   Min.   :  7   Min.   :  3.00   0:13173       N:7457  
##  1st Qu.: 0.8000   1st Qu.: 19   1st Qu.: 17.00   1: 3535       Y:9251  
##  Median : 0.9000   Median : 23   Median : 26.00                         
##  Mean   : 0.8889   Mean   : 26   Mean   : 39.77                         
##  3rd Qu.: 1.0000   3rd Qu.: 28   3rd Qu.: 44.00                         
##  Max.   :10.3000   Max.   :778   Max.   :933.00

#Applicazione trasformate logaritmiche 
test$fasting.blood.sugar<-log(test$fasting.blood.sugar)
test$triglyceride<-log(test$triglyceride)
test$HDL<-log(test$HDL)
test$LDL<-log(test$LDL)
test$serum.creatinine<-log(test$serum.creatinine)
test$AST<-log(test$AST)
test$Gtp<-log(test$Gtp)

#Standardizzazione
test.numeric <- test[,-c(1,5,6,18,19,20)]
for (i in 1:nrow(matrix_indicators)){
  test.numeric[, i] <- (test.numeric[, i] - matrix_indicators[i, "media"])/(matrix_indicators[i, "dv"])
}
test <- cbind(test.numeric, test$gender, test$hearing.left., test$hearing.right.,
               test$dental.caries, test$tartar, test$smoking)
colnames(test)[15:20] <- c('gender','hearing.left.','hearing.right.','dental.caries','tartar','smoking')

#Unione sub.train e validation
train.finale <- rbind(sub.train,validation)

Si riallena il modello sul train completo e si valuta sul test.

set.seed(123)
rf.finale <- randomForest(smoking ~., data=train.finale, cutoff=c(0.7,0.3)) 
rf.pred.class.finale <- predict(rf.finale, test[,-20])
conf5 <- confusionMatrix(rf.pred.class.finale, test$smoking, positive='1', mode='everything')
conf5

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0 7311  269
##          1 3222 5906
##                                           
##                Accuracy : 0.7911          
##                  95% CI : (0.7848, 0.7972)
##     No Information Rate : 0.6304          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.592           
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
##                                           
##             Sensitivity : 0.9564          
##             Specificity : 0.6941          
##          Pos Pred Value : 0.6470          
##          Neg Pred Value : 0.9645          
##               Precision : 0.6470          
##                  Recall : 0.9564          
##                      F1 : 0.7719          
##              Prevalence : 0.3696          
##          Detection Rate : 0.3535          
##    Detection Prevalence : 0.5463          
##       Balanced Accuracy : 0.8253          
##                                           
##        'Positive' Class : 1               
## 

La sensitivity è superiore al 95% e l’F1 score è del 77%, evidenziando un’ottima capacità previsiva da parte del modello scelto per la classe di riferimento. L’accuracy è comunque vicina all’80%, evidenziando una capacità del modello di classificare correttamente su entrambe le classi.

Si grafichi la classificazione con il modello scelto evidenziando le osservazioni misclassificate.

misclassificatetest<-(classError(rf.pred.class.finale, test$smoking)$misclassified)

plotv<-ggplot(test, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point()+
  scale_colour_manual(values=c("darkgreen", "red"))+
  labs(title="Classificazione originale vere label nel test set")+
  theme_bw()
plotv

plotc<-ggplot(test, aes(x=Gtp, y=Cholesterol, col=smoking))+
  geom_point()+
  scale_colour_manual(values=c("darkgreen", "red"), name="classificazione")+
  geom_point(data=test[misclassificatetest,],  col="blue", alpha=0.5)+
  labs(title="Classificazione secondo il metodo random forest nel test set")+
  theme_bw()
plotc

I problemi della classificazione rimangono nella parte centrale. A causa dell’elevato numero di osservazioni la rappresentazione è poco informativa: le unità misclassificate sono in tutto 3491 su un totale di 16708. Inoltre, al fine del problema in analisi la percentuale di misclassificate nella classe di fumatori è ancora più bassa.

7 Conclusioni

Questa analisi si conclude quindi con la scelta di una modellistica che possa essere utile alla previsione dello stato di fumatore (o non fumatore) in un caso assicurativo, in cui i costi associati ai diversi tipi di errore incidono in modo rilevante.

Fonti:

• https://www.kaggle.com/datasets/kukuroo3/body-signal-of-smoking
• https://www.data.go.kr/
• https://www.nhis.or.kr/
• https://www.koreatimes.co.kr/
• https://www.my-personaltrainer.it/
• https://www.epicentro.iss.it/